/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of LINAL (LINear ALgebra)
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <cublas.h>

#define M 6
#define N 5

#define IDX2C(i,j,ld) (((j)*(ld))+(i))

void modify (double *m, int ldm, int n, int p, int q, double alpha,double beta)
{
  cublasDscal (n-p, alpha, &m[IDX2C(p,q,ldm)], ldm);
  cublasDscal (ldm-p, beta, &m[IDX2C(p,q,ldm)], 1);
}

int main(int argc, char *argv[])
{
  int i, j;
  cublasStatus stat;
  double* devPtrA;
  double* a = 0;
  a = (double *)malloc (M * N * sizeof (*a));
  if (!a) {
    printf ("host memory allocation failed");
    return 1;
  }
  for (j = 0; j < N; j++) {
    for (i = 0; i < M; i++) {
      a[IDX2C(i,j,M)] = i * M + j + 1;
    }
  }

  cublasInit();
  stat = cublasAlloc (M*N, sizeof(*a), (void**)&devPtrA);
  if (stat != CUBLAS_STATUS_SUCCESS) {
    printf ("device memory allocation failed");
    return 1;
  }
  cublasSetMatrix (M, N, sizeof(*a), a, M, devPtrA, M);
  modify (devPtrA, M, N, 1, 2, 16.0f, 12.0f);

  for (j = 0; j < N; j++) {
    for (i = 0; i < M; i++) {
      printf ("%8.2f", a[IDX2C(i,j,M)]);
    }
    printf ("\n");
  }
  printf ("\n");

  cublasGetMatrix (M, N, sizeof(*a), devPtrA, M, a, M);


  for (j = 0; j < N; j++) {
    for (i = 0; i < M; i++) {
      printf ("%8.2f", a[IDX2C(i,j,M)]);
    }
    printf ("\n");
  }
  printf ("\n");

  cublasFree (devPtrA);
  cublasShutdown();
  for (j = 0; j < N; j++) {
    for (i = 0; i < M; i++) {
      printf ("%8.2f", a[IDX2C(i,j,M)]);
    }
    printf ("\n");
  }

  return 0;
}
